Dummyvariabler (dummy variables / indicator variables)

• Vi har hittills bara använt numeriska förklaringsvariabler i våra regressionsmodeller

• Det är dock också möjligt att använda kategoriska förklaringsvariabler för att göra gruppvisa regressionsanalyser

• Vi skulle t.ex. kunna dela upp analysen av bränsleförbrukning efter bilens växellåda (manuell/automatisk)

Dummyvariabler (dummy variables / indicator variables)

• Antag att vi vill prediktera amerikanska studenters snittbetyg (GPA) med SAT score som förklaringsvariabel (SAT är ett “amerikanskt högskoleprov”)
• Diagrammet tyder på ett positivt samband mellan dessa numeriska variabler

Motivering av regression med dummyvariabler

• Antag att vi också har en variabel attendance, som har värdet yes om studenten deltagit på minst 75% av föreläsningarna och no i annat fall
• I denna figur är punkter som representerar studenter med hög närvaro gula
• Studenter med hög närvaro tycks generellt ha högre snittbetyg, givet sin SAT-poäng – vi vill gärna använda detta i vår regressionsmodell!

Skapa en dummyvariabel

• Vi har konstaterat att variabeln attendance är kategorisk
• Vi kan använda den i en regressionsmodell om vi kodar om den så att den blir en numerisk dummyvariabel

Skapa en dummyvariabel

GPA$high_attendance_dummy <-
  ifelse(test = (GPA$high_attendance == "yes"), yes = 1, no = 0)
GPA[c(1, 16), ]

    GPA  SAT high_attendance high_attendance_dummy
1  2.40 1714              no                     0
16 3.17 1872             yes                     1

Regression med en dummyvariabel

• Nu när vi har skapat en dummyvariabel använder vi den på samma sätt som en vanlig numerisk variabel
• Vår modell för att prediktera snittbetyg blir

\[\widehat{\text{GPA}} = b_0 + b_1 \cdot \text{SAT} + b_2 \cdot \text{HighAttendance}\]

lm(GPA ~ SAT + high_attendance_dummy, data=GPA)$coefficients

          (Intercept)                   SAT 
          0.643850459           0.001399802 
high_attendance_dummy 
          0.222644088 

Tolkning av regression med dummyvariabel

• Nu när vi har våra koefficienter kan vi skriva modellen som

\[\widehat{\text{GPA}} = 0.644 + 0.0014 \cdot \text{SAT} + 0.223 \cdot \text{HighAttendance}\]

• Tolkning: Modellen uppskattar att en student med hög närvaro har ett snittbetyg som är 0.223 högre än en student som inte har hög närvaro, givet en viss SAT-poäng

Tolkning av regression med dummyvariabel

• Koefficienten som är kopplad till dummyvariabeln visar hur kategorin som kodats till 1 förhåller sig till kategorin som kodats till 0

• När vi kodade hög närvaro som 1 och ej hög närvaro som 0 fick vi \[\widehat{\text{GPA}} = 0.644 + 0.0014 \cdot \text{SAT} + 0.223 \cdot \text{HighAttendance}\]

• Kategori 1 har högre snittpoäng än kategori 0

Tolkning av regression med dummyvariabel

\[\widehat{\text{GPA}} = 0.867 + 0.0014 \cdot \text{SAT} - 0.223 \cdot \text{NotHighAttendance}\]

• Denna modell uppskattar att en student som inte har hög närvaro har ett snittbetyg som är 0.223 lägre än en student som har hög närvaro, givet en viss SAT-poäng

Regression med dummyvariabel

Modellval

• Hittills har vi framför allt undersökt
  • Hur vi tolkar en regressionsmodell
  • Hur vi transformerar variabler för hitta samband som är någorlunda linjära
  • Hur vi använder måttet \(R^2\) som säger oss något om hur väl modellen förklarar responsvariabelns värden

Val av variabler

• En naturlig fråga att ställa är: Varför välja ut vissa variabler och välja bort andra? Kan vi inte bara inkludera alla tillgängliga variabler?
• När vi tar med alla variabler gör vi \(R^2\) så stort som möjligt, men vi vet att \(R^2\) alltid ökar när vi lägger till fler variabler (och inte alltid är att lita på)

Flexibilitet

• Här ser vi 3 regressionsmodeller som estimerar en bils bränsleförbrukning med förklaringsvariabeln hästkrafter
• Vi bryr oss i det här sammanhanget inte om hur modellerna är konstruerade, utan vi konstaterar bara att de har olika stor flexibilitet
• En mer flexibel en modell följer responsvariabelns värden närmare

Flexibilitet

• Den minst flexibla modellen är en rät linje
• Den något mer flexibla modellen i mitten följer det övergripande mönstret
• Den mest flexibla modellen anpassar sig tydligt efter de enskilda datapunkterna

Flexibilitet

Generaliserbarhet

• Vi säger att en modell är generaliserbar om den hade sett ungefär likadan ut om vi anpassade den till andra (men liknande) data
• Vi kan tänka att vi helt plötsligt får data på 32 andra bilmodeller, och skattar de båda modellerna igen
• Vad skulle hända då?

Generaliserbarhet

• Den vänstra modellen skulle förmodligen bli ungefär likadan, då den fångar upp ett övergripande mönster (mer generaliserbar)
• Den högra modellen skulle troligtvis se helt annorlunda ut, då den är väldigt exakt anpassad till de enskilda punkterna (mindre generaliserbar)

Generaliserbarhet och överanpassning

• En flexibel modell med dålig generaliserbarhet sägs vara överanpassad (overfitted)
• Den högra modellen är alltså en överanpassad modell, den är för anpassad till de enskilda punkterna, och bortser så från det övergripande mönstret

Generaliserbarhet och underanpassning

• En modell sägs vara underanpassad (underfitted) om den
  • Inte är flexibel, och
  • Inte fångar det övergripande mönstret i data
• Den röda linjen i bilden beskriver en underpassad modell, som inte fångar det uppenbara bågformade mönstret i datamaterialet

Trade-off mellan över- och underanpassning

• Vi har sett att både över- och underanpassade modeller är problematiska, och målet när vi väljer modell är att hitta en bra medelväg mellan de båda
• Vi vill ta fram en modell som är
  • Tillräckligt flexibel för att fånga det övergripande mönstret
  • Tillräckligt restriktiv för attinte fånga alla slumpvisa variationer

Modellval med hjälp av \(R_{\text{adj}}^2\)

• Vi har tidigare nämnt \(R^2\), och att \(R^2\) alltid växer när vi lägger till variabler
• \(R^2 = 1\) när modellen träffar alla punkter, alltså när den är helt överanpassad
• Om vårt mål är att maximera \(R^2\) kommer vi så att överanpassa modellen

Träningsdata och testdata

• Vi har sagt att en generaliserbar modell är en modell som fungerar även för nya data (data som vi inte använde när vi anpassade modellen)
• Ett bra sätt att bedöma om modellen är generaliserbar är att helt enkelt testa hur bra prediktionerna blir på nya data

Träningsdata och testdata

• Nedan ser vi en model som inte är flexibel, och en modell som är det
• Den flexibla modellen passar träningsdata bra, men nya data väldigt dåligt

Träningsdata och testdata

• Vi kan t.ex. använda 70% av observationerna som träningsdata, och de övriga 30% av observationerna som testdata
• Om observationerna ligger i någon särskild ordning bör vi sortera dem slumpmässigt innan vi gör uppdelningen
• Vi vill att både tränings- och testdata ska vara representativa för hela datamaterialet (t.ex. inte bara ha hushåll från en stadsdel i testdata)

Träningsdata och testdata

• För att utvärdera en modell på testdata följer vi stegen nedan

  1. Dela upp observationerna i träningsdata och testdata
  2. Anpassa regressionsmodellen med hjälp av träningsdata
  3. Utvärdera modellen med hjälp av testdata

Modellval - träningsdata och testdata

• Vi räknar ut RMSE på testdata med formeln \[\text{RMSE}_\text{test} = \sqrt{\cfrac{\sum_{i=1}^{n_\text{test}} (y_i - \hat{y}_i)^2}{n_\text{test}}},\]

• Vi summerar över alla observationer i testdata, och \(n_\text{test}\) är så antalet observationer i testdata

Träningsdata och testdata

• Vi kommer ihåg att SSE (Sum of Squared Errors) betecknar summan av kvadrerade prediktionsfel, \[\text{RMSE}_\text{test} = \sum_{i=1}^{n_\text{test}} (y_i - \hat{y}_i)^2\]

• Vi kan alltså skriva RMSE som \[\text{RMSE}_\text{test} = \sqrt{\cfrac{SSE_\text{test}}{n_\text{test}}}\]

Modellval - träningsdata och testdata

• Det är ofta lättast (rent räknemässigt) att beräkna \(\text{RMSE}_\text{test}\) i flera steg
  1. Beräkna \(\text{SSE}_\text{test} = \sum_{i=1}^{n_\text{test}} (y_i - \hat{y}_i)^2\)
  2. Beräkna \(\text{MSE}_\text{test} = \cfrac{\text{SSE}_\text{test}}{n_\text{test}}\)
  3. Beräkna \(\text{RMSE}_\text{test} = \sqrt{\text{MSE}_\text{test}}\)
• Vi använder ofta \(\text{SSE}_\text{test}\) till flera olika saker, så det är bra att beräkna det separat för att slippa repetera

Träningsdata och testdata

• Vi tittar nu på ett exempel, där vi
  1. Delar upp ett dataset i träningsdata och testdata
  2. Anpassar modellen med hjälp av vår träningsdata
  3. Utvärderar modellen genom att räkna ut RMSE för vår testdata

Exempel

• Vi använder ett dataset med variablerna \(y\), \(x_1\), \(x_2\) och \(x_3\) – där \(y\) är responsvariabel och övriga variabler är förklaringsvariabler
• Vilka förklaringsvariabler vi inkluderar är upp till oss själva, och målet är att hitta en modell som ger bra prediktioner i vår testdata
• Vi har totalt 100 observationer, och de första ser ut såhär:

         y       x1       x2         x3
1 31.41966 16.88882 5.599934 -0.7104066
2 32.39954 25.40119 3.996941  0.2568837
3 26.59989 18.95261 4.931678 -0.2466919
4 39.44798 27.01130 7.726843 -0.3475426

Steg 1: Dela upp data

• Vi börjar med att dela in våra data i en träningsdel och ett testdel
• Det gör vi genom att först lägga observationerna i slumpvis ordning
• Därefter låter vi de första 70 observationerna utgöra träningsdata
• Till sist får de återstående 30 observationerna bli testdata

library(dplyr)
data_randomorder <- data %>% slice_sample(prop=1)
datatrain <- data_randomorder %>% slice(1:70)
datatest <- data_randomorder %>% slice(71:100)

Steg 2: Träna modellen \(y = b_0 + b_1 x_1 + b_2 x_2\)

model1 <- lm(y ~ x1 + x2, data=datatrain) #Träna modellen
summary(model1) #Se sammanfattande resultat

Steg 3: Utvärdera på testdata

• Vi vill nu beräkna \(\text{RMSE}_\text{test}\), och beräknar \(\text{SSE}_\text{test}\), \(\text{MSE}_\text{test}\) och \(\text{RMSE}_\text{test}\) i tur och ordning

rsquared1 <- summary(model1)$r.squared
y_hatt_test <- predict(model1, newdata=datatest)
y_test <- datatest$y
SSE <- sum((y_test - y_hatt_test)^2)
n_test <- 30
MSE <- SSE / n_test
RMSE <- sqrt(MSE)
sprintf("R-squared=%.7f, SSE=%.7f", rsquared1, SSE)

[1] "R-squared=0.8709888, SSE=146.0526235"

sprintf("MSE=%.7f, RMSE=%.7f", MSE, RMSE)

[1] "MSE=4.8684208, RMSE=2.2064498"

Upprepa för modellen \(y = b_0 + b_1 x_1 + b_2 x_2 + b_3 x_3\)

model2 <- lm(y ~ x1 + x2 + x3, data=datatrain)
rsquared2 <- summary(model2)$r.squared
y_hatt_test <- predict(model2, newdata=datatest)
y_test <- datatest$y
SSE2 <- sum((y_test - y_hatt_test)^2)
n_test <- 30
MSE2 <- SSE2 / n_test
RMSE2 <- sqrt(MSE2)
sprintf("R-squared=%.7f, SSE=%.7f", rsquared2, SSE2)

[1] "R-squared=0.8710176, SSE=146.5228725"

sprintf("MSE=%.7f, RMSE=%.7f", MSE2, RMSE2)

[1] "MSE=4.8840958, RMSE=2.2099990"

Jämförelse av modellerna

• Vi jämför modellerna i tabellen nedan, där vi har
  • Modell 1: \(y = b_0 + b_1 x_1 + b_2 x_2\)
  • Modell 2: \(y = b_0 + b_1 x_1 + b_2 x_2 + b_3 x_3\)

Träningsdata och testdata

• Att dela upp ett dataset i träningsdata och testdata ger oss en bättre bild av hur användbar en modell är
• Vi vet dock inte säkert om bilden blir helt rättvis, särskilt om det är ett dataset med få observationer

Träningsdata och testdata

• Vi kodar upp exemplet i R, och börjar med de 12 sista observationerna som testdata

# De 12 sista observationerna används i vårt testset
carstrain <- mtcars %>% slice(1:20)
carstest <- mtcars %>% slice(21:32)
lmod1 <- lm(litermil ~ viktton, data=carstrain)
y_hatt <- predict(lmod1, newdata=carstest)
RMSE <- sqrt(mean((carstest$litermil - y_hatt)^2))
print(RMSE)

[1] 0.1984769

Träningsdata och testdata

• Vi gör nu om analysen, fast med de 12 första observationerna som testdata

# De 12 första observationerna används i vårt testset
carstrain <- mtcars %>% slice(13:32)
carstest <- mtcars %>% slice(1:12)
lmod1 <- lm(litermil ~ viktton, data=carstrain)
y_hatt <- predict(lmod1, newdata=carstest)
RMSE <- sqrt(mean((carstest$litermil - y_hatt)^2))
print(RMSE)

[1] 0.1790871

Träningsdata och testdata

• Med de sista 12 observationerna som testdata fick vi \(\text{RMSE}=0.198\)
• Med de första 12 observationerna som testdata fick vi \(\text{RMSE}=0.179\)
• Vi ser så att \(\text{RMSE}\) skiljer sig relativt mycket, beroende på hur vi delar upp våra data
• Ett sätt att komma runt detta problem är att använda något som kallas för korsvalidering (cross validation)

Korsvalidering

• Korsvalidering är en metod som låter oss använda alla observationer både som träningsdata och som testdata
• När vi använder korsvalidering delar vi upp våra data i tränings- respektive testdelar flera gånger
• Efter varje uppdelningen får vi en ny uppsättning punkter som testdata
• Varje observation hamnar i testdata en gång, och i träningsdata i övriga uppdelningar

Korsvalidering

• Varje uppdelning kallas för en fold
• Korsvalisering kallas ibland \(K\)-fold cross validation, där \(K\) är antalet folds
• Korsvalidering med \(K = 4\) betyder exempelvis att vi delar in vår data i tränings- och testdata på 4 olika sätt

Korsvalidering

• För var och en av de 4 uppdelningarna räknar vi ut SSE på testdata
• Vi använder sedan summan av de 4 SSE-värdena för att beräkna RMSE

Korsvalidering

• Antag att vi har ett dataset med 100 observationer och \(K=4\)
  • För \(k=1\) gäller då att observation 1 till 25 utgör testdata
  • För \(k=2\) gäller att observation 26 till 50 utgör testdata,
  • Och så vidare

Korsvalidering

• Vårt mål är att räkna ut RMSE med hjälp av våra 4 folds
• Vi börjar med att räkna ut SSE för varje enskild fold, och summerar sedan våra fyra värden som \[\text{SSE}_\text{cv} = \text{SSE}_\text{test}^{(1)} + \text{SSE}_\text{test}^{(2)} + \text{SSE}_\text{test}^{(3)} + \text{SSE}_\text{test}^{(4)}\]

Korsvalidering

• Procedur för enkel uppdelning i tränings- och testdata:
  1. Räkna ut \(\text{SSE}_\text{test}\)
  2. Räkna ut \(\text{MSE}_\text{test} = \text{SSE}_\text{test} / n_{\text{test}}\)
  3. Räkna ut \(\text{RMSE}_\text{test} = \sqrt{\text{MSE}_\text{test}}\)

Exempel på korsvalidering med \(K=2\)

• Vi arbetar nu igenom ett exempel på korsvalidering, med hjälp av R
• Vi använder våra bildata med 32 observationer, och utvärderar den enkla modellen med vikt som förklaringsvariabel
• Vi sätter \(K=2\), det vill säga vi använder 2 folds

Exempel på korsvalidering med \(K=2\)

• Med koden nedan kan vi beräkna \(\text{SSE}_\text{test}^{(1)}\) och \(\text{SSE}_\text{test}^{(2)}\)

# Räkna ut SSE(1)
carstest <- mtcars %>% slice(1:16)
carstrain <- mtcars %>% slice(17:32)
lmod1 <- lm(litermil ~ viktton, data=carstrain)
y_hatt1 <- predict(lmod1, newdata=carstest)
y1 <- carstest$litermil
SSE1 <- sum((y1 - y_hatt1)^2)

# Räkna ut SSE(2)
carstest <- mtcars %>% slice(17:32)
carstrain <- mtcars %>% slice(1:16)
lmod2 <- lm(litermil ~ viktton, data=carstrain)
y_hatt2 <- predict(lmod2, newdata=carstest)
y2 <- carstest$litermil
SSE2 <- sum((y2 - y_hatt2)^2)

Exempel på korsvalidering med \(K=2\)

Vi skriver ut värdena på \(\text{SSE}_\text{test}^{(1)}\) och \(\text{SSE}_\text{test}^{(2)}\)

#Skriv ut värdet på SSE(1) och SSE(2)
sprintf("SSE1 = %.4f, SSE2 = %.4f", SSE1, SSE2)

[1] "SSE1 = 0.5331, SSE2 = 1.0270"

Exempel på korsvalidering med \(K=2\)

• Bilderna visar våra data och de två regressionsmodellerna (en för varje fold)
• Om vi kvadrerar de röda linjerna i bilderna till höger och summerar kvadraterna får vi \(\text{SSE}_\text{test}^{(1)}\) (övre) och \(\text{SSE}_\text{test}^{(2)}\) (undre)

Exempel på korsvalidering med \(K=2\)

• Vi har räknat ut
  • \(\text{SSE}_\text{test}^{(1)} = 0.5331\)
  • \(\text{SSE}_\text{test}^{(2)} = 1.0270\)
• Nu kan vi räkna ut

\[\text{SSE}_\text{CV} = \text{SSE}_\text{test}^{(1)} + \text{SSE}_\text{test}^{(2)} = 0.5331 + 1.0270 = 1.5601\]

\[\text{MSE}_\text{CV} = \cfrac{\text{SSE}_\text{CV}}{n} = \cfrac{1.5601}{32} = 0.04875312\] \[\text{RMSE}_\text{CV} = \sqrt{\text{MSE}_\text{CV}} = \sqrt{0.04875312}=0.2208011\]

Exempel på korsvalidering med \(K=2\)

• Normalt gör vi våra uträkningar i programkoden

n <- 32
SSE_cv <- SSE1 + SSE2
MSE_cv <- SSE_cv / n
RMSE_cv <- sqrt(MSE_cv)
sprintf("RMSE är %.5f", RMSE_cv)

[1] "RMSE är 0.22080"

Exempel på korsvalidering med \(K=2\)

• Nu har vi räknat ut \(\text{RMSE}_\text{CV}\) för en regressionsmodell
• För att jämföra flera modeller måste vi räkna ut \(\text{RMSE}_\text{CV}\) för varje modell
• När vi har gjort detta väljer vi den modell om har lägst \(\text{RMSE}_\text{CV}\)

Credits

Dessa slides skapades av Karl Sigfrid för kursen Statistik och Dataanalys I och har uppdaterats av Oskar Gustafsson och Valentin Zulj